% 第一题
clear;
location=readmatrix('附件');% location为定日镜xy坐标
% figure;plot(location(:,1),location(:,2),'.');
% title('附件定日镜坐标信息');axis equal;
loc_jire=[0,0,80];%集热器中心坐标
h=4;% 定日镜高度
loc_dingri=[location,h+zeros(size(location,1),1)]; %定日器中心坐标
s_reflect=loc_jire-loc_dingri; % 反射光线的方向向量
s_reflect=s_reflect./sqrt(s_reflect(:,1).^2+s_reflect(:,2).^2+s_reflect(:,3).^2);%单位化
D = [306;337;0;31;61;92;122;153;184;214;245;275];
% D为以春分作为第0天起算的天数 12个月
ST=9:1.5:15;% 当地时间
% i=D长12 j=ST长5 共60次
alphas=zeros(12,5);gamas=zeros(12,5);s_in=zeros(12,15); %初始化
phi=39.4*pi/180;% 当地纬度
for i=1:12 %第i行为日期
for j=1:3 %第j列为时间
delta(i,j)=asin( sin(2*pi*D(i)/365)*sin(2*pi*23.45/360) );% 太阳赤纬角
w(i,j)=(ST(j)-12)*pi/12;% 太阳时角
alphas(i,j)=( asin( cos(delta(i,j))*cos(phi)*cos(w(i,j))+sin(delta(i,j))*sin(phi) ));% 太阳高度角
gamas(i,j)=real(acos((sin(delta(i,j))-sin(alphas(i,j)))*sin(phi))/(cos(alphas(i,j))*cos(phi)));
%太阳方位角
% if gamas(i,j)>pi
% gamas(i,j)=2*pi-gamas(i,j);
% end
s_in(i,3*(j-1)+(1:3))=-[sin(gamas(i,j)),cos(gamas(i,j)),tan(alphas(i,j))];%
s_in(i,3*(j-1)+(1:3))=s_in(i,3*(j-1)+(1:3))./norm(s_in(i,3*(j-1)+(1:3)));% 单位化
end
end
%fprintf('%f\n', (gamas));

s_in(:,10)=-s_in(:,4);
s_in(:,11)=s_in(:,5);
s_in(:,12)=s_in(:,6);
s_in(:,13)=-s_in(:,1);
s_in(:,14)=s_in(:,2);
s_in(:,15)=s_in(:,3);
% i=D长12 j=ST长5 共60次
n_dingri=zeros(size(location,1),3); % 第i天,第j时,定日器法向量
%% 余弦效率
eta_cos=zeros(12,5,size(location,1));
for i=1:12 %第i行为日期
for j=1:5 %第j列为时间
n_dingri=s_in(i,3*(j-1)+(1:3))-s_reflect;% 确定每个定日器的法向量
n_dingri=n_dingri./sqrt(n_dingri(:,1).^2+n_dingri(:,2).^2+n_dingri(:,3).^2);
for k=1:size(location,1)
eta_cos(i,j,k)=abs(dot(n_dingri(k,:),s_in(i,3*(j-1)+(1:3))));
end
end
end
month_eta_cos=zeros(12,1);
for i=1:12
month_eta_cos(i)=sum(eta_cos(i,:,:),'all')/(5*size(location,1));
end
year_eta_cos=mean(month_eta_cos);
fprintf('结果：%f',year_eta_cos);